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Abstract 

In the postgenome era many efforts have been dedicated to systematically elucidate the complex 
web of interacting genes and proteins. These efforts include experimental and computational 
methods. Microarray technology offers an opportunity for monitoring gene expression level at the 
genome scale. By recourse to information theory, this study proposes a mathematical approach to 
reconstruct gene regulatory networks at coarse-grain level from high throughput gene expression 
data. The method provides the a posteriori probability that a given gene regulates positively, 
negatively or does not regulate each one of the network genes. This approach also allows the 
introduction of prior knowledge and the quantification of the information gain from experimental 
data used in the inference procedure. This information gain can be used to chose genes to be 
perturbed in subsequent experiments in order to refine the knowledge about the architecture of an 
underlying gene regulatory network. The performance of the proposed approach has been studied 
by in numero experiments. Our results suggest that the approach is suitable for focusing on size- 
limited problems, such as, recovering a small subnetwork of interest by performing perturbation 
over selected genes. 
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Introduction 



Gene expression is regulated by proteins that enhance or block polymerase binding at 
the promoter region. These biochemical reactions constitute the edges of the gene regula- 
tory networks. One of the key issues in modern biology is the elucidation of the structure 
and function of gene regulatory circuits at the system level To address this challenge 
many efforts have been devoted to the task of developing computational methods capable of 
inferring the interaction between genes from expression levels both on small pathways 
as on genome- wide scale (see J4] for a review). Several models for gene regulatory networks 
rave been proposed in order to infer network interactions la, such as Bayesian networks 

inn n nnn 

3, |8|, |9j, Boolean networks [10(, linear model [Uj, [12|, [13] . Once a regulatory network model 



has been chose, it is possible, in principle, to recover its parameters with some accuracy. Of 
course, more detailed models will require more extensive experimental data. In general this 
data is not available for the genome-wide scale assuming complex model. However, we can 
concentrate on simpler task, such as: who is regulating whom? and, Is that an up-regulation 
or a down-regulation? The idea behind restricting our questions at this qualitative infor- 
mation level, is to reduce the amount of data needed to infer valuable and robust biological 
knowledge even when dealing with noise data. In any case, the detailed information offered 
by more detailed modeling is not useful without a careful significance analysis of these pre- 
dictions. In this sense, this study proposes a mathematical approach to infer gene networks 
at the coarse grain level. The inference process is to be accomplished according to Ockham's 
razor, i.e., with the minimum number of assumptions compatible with the available data. To 
do that, the information theory (IT) is used within the framework of the maximum entropy 



principle 14], [ijj] . IT has proved to be of uti 



expression and network reconstruction 16|, [l7| , where gene expression levels were regarded 



ity in devising techniques for analyzing gene 



as random variables. Here, complementing these previous works, each putative interaction 
has been considered as a random variable. In numero experiments show that, in this case, 
the IT parlance also provides a powerful framework to discuss questions related to the mod- 
eling process such as: (i) how to incorporate a priori information about the gene interaction; 
(ii) how to asses the likelihood of the inferred paths; (iii) how to quantify the information 
provided by the experimental data; and (iv) how to design experiments in order to identify 
subnetworks. 
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The IT approach 



In general a genetic network can be modeled by a set of non-linear differential equations 
x % — fi ( x i(t)i ■ ■ ■ 5 x N{t)), where Xi(t) is the expression level of gene i at time t and / f is the 



regulatory function governing the expression of gene % 18j . Near a steady state the nonlinear 



system can be approximated by a set of linear differential equations, x = Wx where W, is 



a weighted connectivity matrix 19J. In order to uncover the connectivity matrix, we can 
apply a stimulus b = (fei(t), . . . , fejv(t)) T , then measure simultaneously the mRNA levels 
relative to N genes x and estimate the derivative x. Repeating the procedure M times 
we get a measurement matrix X where columns denote the experiments and where rows 
indicate individual genes. Thus, we can approximate the dynamics by 

X = WX + B (1) 

where X and B follow the same notation as X. 

Usually, inferring genetic network attempts to retrieve the weight matrix W, where the 
elements describe the type and strength of influence of gene j on gene i (wij > indi- 
cates activation, > indicates repression, and a zero indicates no influence). However, 
without a careful assessment of the significance of the weight-values, this could lead to the 
conclusion that the network is fully connected in contradiction to the well-known fact that 
gene regulatory networks are sparse networks. 

In the present work, the maximum entropy principle is applied to obtain the probability 
distribution from the data Dm = {x, X, B j, over the possible matrix W. After that, using 
maximum likelihood criterion, the gene interaction matrix I is selected. The elements 7y 
can take only three values, depending on the type of influence of gene j on gene i, Iij = 1 
for activation (direct or indirect), 1^ = —1 for repression and 1^ = when gene j does not 
have influence on gene i. In order to infer weights consistent with Dm, it is assumed that 
each set of weights W is realized with probability P (W\Dm)- In other words, a normalized 
probability distribution is introduced over the possible sets W, which satisfy 

(W) = J P(W\D M )WdW. (2) 

The relative entropy related to an a priori probability distribution P , is given by 

~P(W\D M ) 



H r (D M \P ) = -Jp(W\D M ) In 



P (W) 



dW, (3) 
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where Pq (W) is an appropriate a priori distribution. The negative relative entropy H r , 



known as Kullback-Leibler distance 



20|, defines the information gained after Dm has been 



used in the inference procedure. Thus, in this framework, the inference process takes place 
through a modification of the probability distribution on weights space due to incoming 
data. 

Among all possible distributions P(W\Dm) consistent with Dm, P (W\Dm) has been 
selected which comprises no unjustified prejudice. Thus, following the central tenets of the 
maximum entropy principle, relative entropy is maximized subject to the constraints Eq. El 
Thus, the a posteriori probability distribution yields, 



P (W\D M ) = exp (- (1 + A )) exp (-W ■ T) P (W) 



(4) 



where Ao is Lagrange multiplier associated to the normalization condition, and T the La- 
grange multipliers associated to the constraints Eq. El which are determined once Pq is 
properly selected. 

In order to select Pq, it is assumed that the weights are restricted to the values of i.e. 



W; 



0, ±1 and then a three-peaked a priori distribution is used, which is described by 



N 



P (W) 



(^ J + 1 ) 2 



Pij e 2a +P?i e 



' + 



(5) 



where pf, is the a priori probability for gene j to regulate positively (x = +), negatively 
(x = — ) or to not regulate (x = 0) gene i. Of course + pfj + Pij = 1 for each pair i,j. The 
parameter a can be regarded as a constraint smoothness parameter. Replacing this choice 
in Eq. H] the a posteriori probability distribution is obtained as a sum of three Gaussians, 



P(W\D 



1 



Mi 



+Pij e ~ 



(2ira) N /'< 



N 

n 



pi e_ 



+ Pij e 



(6) 



where pf, is the a posteriori probability for gene j regulate positively (x = +), negatively (x = 
— ) or to not regulate (x = 0) gene i. These probabilities are defined by p% = p^e _r ^ / , 
Pij = Pij eVij l z ii an d Pij — Pij/ Z ij, where z {j = 1 + p J (e~ riJ - + p^ (e r ^ - guarantee 
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|r| - In ( Zij ) 



normalization. Furthermore, the relative entropy of the a posteriori distribution Eq. [3] is 
given by 

N 

H r (D M , P ) = - E h {i\D M , Po) , (7) 

j 

where is the information gain of gene i with respect to Po obtained from using the data 
D M which is defined by 

N 

I g (i\D M ,P ) = E 

3 

— O'/; 1 ' ! - /'-,!'-,' '' ') • (8) 

Zij 

The multipliers 1^- are obtained after solving the equation 

(wu) = al„ + (p+ t 1 - p".e r -) . (9) 

where (■?%) are subject to the constraints imposed by Dm- Our central idea is that of 
reinterpreting, following the information in Dm in a particular fashion, 

X-B = (W)X. (10) 

Thus, all of the possible networks that are consistent with Eq. [TO] can be written as 

(W) = (X - B) ■ U ■ diag(sT 1 ) ■ V T + C ■ V T (11) 



C = (Qj) is an iV x iV matrix, where is zero if Sj ^ and is otherwise an arbitrary scalar 
coefficient. U, S and V correspond to the singular value decomposition of matrix X T , i.e. 
X T = U ■ S ■ V T where U is a unitary M x N matrix of left eigenvectors, S is diagonal 
N x N matrix containing the eigenvalues {si, . . . , sn}, and V is a unitary N x N matrix 
of right eigenvectors. Without loss of generality, let all non-zero elements of Sj be listed at 
the end and sj 1 in Eq. [11] are taken to be zero if Sj = 0. The general solution [11] can be 
written as 

(W) = W L2 + C • V T (12) 

where W^ 2 is the particular solution with the smallest L2 norm. If M < N, many weights 
W are compatible with the available information. The information contained in the data set 
Dm can be used in different ways. Each of these leads to a different probability distribution 
which exhibits diverse properties. In this sense, following the prescription (W) = in Eq. 



the knowledge that gene regulatory networks are sparse can be made use of. Thus, we 
have C • V T = — W/, 2 , which is an overdetermined problem 19]. This particular solution 
will be denoted as W Ll . Of course the T is obtained solving Eq. fusing (W) = W La or 
(W) = Wi r In the following sections these alternatives will be considered independently. 
Notice that for M > N, W La = W Ll . 

After determining the a posteriori distribution, the gene interaction matrix / must be 
selected. In order to do that, the maximum likelihood criterion is taken into account, i.e. 
the selection is accomplished choosing the highest a posteriori probability from {p®j,pfj,Pij} 
for each pair i, j. For example if pfj is greater than p^ and p~j, then 1^ = 1 indicating that 
gene j activates the gene i. 

In order to achieve the best model, the idea is to use the information contained in Dm 
and the knowledge that gene regulatory networks are sparse. The formalism presented here 
offers an alternative to the prescription which selects from all possible solutions [TTJ 
This alternative consists in setting pfj = p^ <C In this way the knowledge that gene 
regulatory network is sparse can be introduced by assigning a much lower value to the a 
priori probabilities of interaction than the a priori probabilities of absence of interaction. 
Furthermore, as the inference processes occur row by row, any other relevant a priori in- 
formation about the gene in consideration (such as known interactions, type of gene, etc.) 
could be included in these probabilities. For example, if gene k encode a helix-turn-helix or 
a zinc finger protein, high probabilities can be assigned for column k (pf k and p^ k ). 



Results 



In order to systematically benchmark the inference performance of this method, a linear 
data-generating model was used. The M random inputs (the columns of matrix X) were 
generated in the range [-1,1] and was computed W • X as the system response, where W 
is the matrix to be reconstructed. Thus pairs X, Y constitute the available information 
Dm- In the simulation, it was observed that the mean performance depends on size and 
the degree of connectivity and not on the network type. For this reason, random sparse 
linear networks will be considered, where each gene has k entries in average. To build the 
connectivity matrix W, following procedure was used: for each matrix element a random 
number r between (0,1) was sorted, if r < k/2N a negative random value chosen from a 
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uniform distribution was assigned to the matrix element, if r > l — k/2N the matrix element 
was a positive random number, and otherwise the matrix element was zero. The condition 
k <^ N ensures sparseness. 

By using singular value decomposition and interior point method for L\ regression 
was computed. Subsequently the set of uncoupled nonlinear equation was solved and the a 
posteriori probability for each putative interaction was evaluated. After this procedure the 
most likelihood I can be selected. The performance of the inference procedure was measured 
by the prediction error e = N~ 2 e ij, where is defined by 



Figure 1 depicts the prediction error e as a function of a defined as the ratio of number 
experiments and number of genes, i.e. a = M/N. These have been tested in three different 
size networks with k/N = 0.05, in which all a priori probabilities are assumed to be equals 
(i.e. pfj = p~j = p®j = 1/3) and a = 0.01. For small values of M the method mistakenly 
infers a percentage of interaction which depends on the network size N and k. However, the 
prediction error decays rapidly as a increases and the gene interaction matrix is completely 
recovered with a a value that decreases with the network size. This performance was ob- 
tained using Wii prescription. Similar simulations (data not shown) performed with the 
W^ 2 prescription, reveal that in these cases the prediction error e remains close to unit until 
a — 1, where they decay abruptly. Dependence of performance on the network topology has 
not been detected and similar results were obtained for scale free networks which have more 
biological appeal than the random networks used here. However, these simulations present 
a greater error bar due to the fact that the network building algorithm used here does not 
make networks with a uniform node degree. 

Many times, when dealing with an incomplete data set M <C N, only a percentage of the 
interactions is inferred correctly. If the likelihood of the inferred paths cannot be assessed, 
this partial reconstruction has small predictive value in real life. The methodology proposed 
here can assess the likelihood of the predicted interaction straightforwardly through the a 
posteriori probability. In this sense, only those predicted interactions with an a posteriori 
probability which is greater than some significance level can be selected. To illustrate this 
issue, a network with 60 genes with k/N = 0.05 was simulated. The related connectivity 
matrix W is represented in Fig. 2 (left), row i corresponds to the genes that regulate the 
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activity of gene i, while column j corresponds to the genes regulated by gene j. The weight 
values Wij are depicted following a linear gray scale, where white (black) corresponds to the 
maximum(minimum) values of weights and the gray background represents the absence of 
interaction. This network is random perturbed in 24 different experiments (a = 0.4). With 
this amount of data usually about ~ 99.5% of the interactions are predicted correctly (see 
Fig. 1). Nevertheless, which interactions were inferred correctly and which were inferred 
wrongly is unknown. By mean of the information theory approach, the a posteriori prob- 
abilities were computed and the inferred interaction matrix I and the associated likelihood 
were derived. Fig. 2 (right) represents the inferred connectivity matrix I, by assuming that 
all a priori probabilities are equal (i.e. pfj = p~j = p\- = 1/3). Red circles indicate wrong 
predictions (1% of the interactions), while green circles indicate the interactions with a a 
posteriori probability greater than 0.99. In this case there are 98(108) interactions where 
the maximum a posteriori probability is greater 0.99(0.95). A more detailed study (data 
not show) revealed that interactions related to higher weight values, are associated to high 
a posteriori probability values. These results suggest that gene networks can be partially 
recovered even with small amounts of data, mainly for those genes that interact strongly. 

Unfortunately all measurements are subject to observational noise, consequently it is 
important to asses to what extent the performance of the inference procedure is affected by 
noise. To simulate this condition in the numerical experiment, the available information Dm 
(both input and output) was corrupted by an additive Gaussian noise with mean zero and 
standard deviation rj. This inference procedure was performed for networks with N = 60, 
in the same condition as for the previous assessment (pfj = p~j = p\j = 1/3 and a = 0.01). 
However, in this case the method based on the prescription of sparseness assumed in Wn 
could not correctly recovery the gene interaction matrix I when the noise level was rj = 0.3 
(even for smaller 77). Figure 3(top) indicates the prediction error by using both W L1 and W L2 
assuming that the a priori probability for activation, repression or absence of interaction 
are equal. This clearly shows that the prediction power decreases as more data becomes 
available. 

However, the network can be partially reconstructed by using an alternative constraint of 
sparsity. This alternative consists in introducing the knowledge of sparseness of the matrix 
through the a priori probabilities. That is achieved by setting pfj m in the inference 
procedure. Fig. 3 (middle) depicts the prediction error as a function of a when the a priori 
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probabilities were set to pfj = 0.025. The sum of these probability values corresponds to 
the percentage of genes that are regulated by one gene. With such a priori information, 
it is possible to reconstruct almost the complete structure of the network (around 95% of 
edges) using more experiments than the number of genes, a ~ 2. The mean node degree 
of the network is generally not known in advance. Notwithstanding, the prediction ability 
is robust for underestimations of the a priori probabilities. Figure 3 (bottom) depicts the 
prediction error as function of a when the a priori probabilities were set to pfj = 0.01. The 
result is almost the same as the previous one. This implies that it is possible to partially 
recover the interaction matrix even with noise data, by setting low values for the a priori 
probabilities pfy In the last two cases, the prediction performance obtained by the ~Wl2 
prescription is comparable with that obtained by W^i using pfj ks prescription, in contrast 
to the case which deals with clean data. Furthermore, when data are corrupted by noise, it 
was observed that prediction error has a peak around a — 1, this peak arises because some 
singular values, associated to the SVD, take small values as consequence of noise. 

The partial recovery referenced above does not pursue recover a closed subnetwork, which 
mainly infer strong interactions around the whole network. However, in many cases this is 
crucial to recover the complete subnetwork associated to a given gene or path of interest. 
The inference approach and information gain tool presented in this study, could be used 
to establish new relationships between genes and to propose new experiments. By means 
of cycles of experiments-datamining, the knowledge about the subnetwork can be refined 
until its complete recovery, even in presence of observational noise. For that purpose the 
following protocol could be used: i) perform an initial perturbation where the gene of interest 
is overexpressed, and obtain the genome expression profile; ii) compute the information gain 
for each gene with this experimental data; iii) select the genes for which the information 
gain is greater than a given threshold; iv) iterate first two steps perturbing each one of the 
genes which were selected in the third step and which have still no been perturbed, until 
no new gene has an information gain greater than the threshold. Figure 4A illustrates the 
result of three of the experiments-datamining cycles. Firstly, the gene which belongs to the 
subnetwork of interest, gene gl, is initially overexpressed (level of 10.0 while the other gene 
levels are random in the range [-0.5,0.5]), then the input-output network is measured, this 
measurement is subject to observational noise with r\ = 0.30. The information gain of this 
experiment is computed for each gene using pf, = = 0.01 as an a priori probability. 
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Subsequently, those genes with I g greater than 1.0 are selected. I g suggests that gene g6 is 
regulated by gl. Repeating the above step with gene g6, the results indicate that genes g2, 
g3 are regulated by g6. The above step is repeated with gene g2 and subsequent genes with 
high information gain values in ensuing experiments, until no new gene with an information 
gain greater than threshold appears. Fig. 4B illustrates a list of experiments where the 
first column corresponds to the gene that was perturbed in the experiment, and the second 
column corresponds to the genes which appear to be regulated by the perturbed gene. In 
the last two experiments no new regulated genes appeared (which were not indicated in the 
first column list). The above analysis provides a causal link between two genes, but it does 
not indicate if the regulation is positive or negative. In order to extract this information, 
the inference analysis was performed using the ten "overexpression experiments" pooled 
in Dm (M = 10). When the inference procedure was applied with this data, 19 out of 
24 interactions in the subnetwork were inferred correctly, 10 of them with an posteriori 
probability greater than 0.99. However, the a priori probabilities provided by the information 
contained in list of Fig. 4B are included, setting pfj = p~j = 0.5 (or 1/3) for all the pairs 
i,j indicated in the list, and pfj = p~j = 0.01 otherwise, 23 out of 24 interactions in the 
subnetwork are inferred, 19 of them with a posteriori probability greater than 0.99, Fig. 4C. 
The performance above obtained does not differ if the inference procedure is implemented 
using W^2 or W^i prescription, of course Wl 2 is computationally cheaper than W^i which 
requires linear programming optimization. 

The above example about subnetwork inference suggests that this novel scheme can be 
re-used regarding further subnetworks until the whole network is recovered with M ~ N 
experiments. 

Discussion and Conclusions 

A novel approach for regulatory network inference is presented in this study. Differently to 
other methods, this approach pursues to infer the type of interaction rather than a weight 
which characterizes the interaction quantitatively. Three main features of the proposed 
method are pointed out. First, it allows to introduce global a priori information about the 
network, as sparseness, and other gene dependent available information, as illustrated in the 
last example Fig. 4C. Second, the information theory formalism provides a way to quantify 
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the likelihood of the inferred paths, by using the a posteriori probabilities computed with the 
method. Last, but not least, information theory formalism also quantifies the information 
gained with the set of data to be used in the inference procedure. 

Furthermore, the IT approach seems to offer promising perspective as a network inference 
protocol; the methodology presented here introduces an information gain measure as a 
bonus. The way in which this quantity could be a useful tool to identify the downstream 
regulated genes in overexpression experiments is illustrated in this study. This feature allows 
a datamining-assisted way of uncovering the whole network with a number of experiments 
equal to the number of genes, even when dealing with a high level of observational noise. 
This IT approach enables the effective use of all the available information, in which each 
experiment is used as an individual constraint. Thus, the ensuing observation level becomes 
much richer than the standard one, where all data define a fitness function to be optimized. 
Efficient management leads to more realistic results in inference. 

The learning protocol presented here constitutes an additional inference technique, which 
should be of interest not only for basic research but also as an application to many interesting 
real world problems without paying an excessive computational cost. 
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FIG. 1: (Color online) Performance. Prediction error e as a function of the ratio a = M/N for 
gene networks with 60 genes (squares), 120 genes (circles) and 240 genes (triangles), averaged over 
50 networks. In all cases the performances were obtained using Wli prescription, equal a priori 



probabilities (i.e. pf^ = p^ = p®j = 1/3 for all i and j), k/N = 0.05 and a = 0.01. 
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FIG. 2: (Color online) Likelihood assessment. Left: connectivity matrix W representation related 
to a random network of 60 genes with k/N = 0.05. Rows correspond to regulated genes, while 
columns correspond to the genes acting as regulators. The interaction weights Wij are represented 
following a linear gray scale, where white corresponds to Wij = 2, while black to Wij = —2. 
Gray background represents the absence of interaction, i.e. Wij = 0. Right: gene interaction 
matrix I inferred after 24 random perturbation experiments, using Wj_,\ prescription, a = 0.01 and 
pfj = Pij = Pij = 1/3- Circles (green in the online figure) indicate the 94 interactions with an a 
posteriori probability greater than 0.99. Wrong predictions (35, ~ 1% of the putative interactions, 
which in this case) correspond to the regulatory inputs of two genes. 
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FIG. 3: (Color online) Inferring with noisy data. Prediction error e as a function of the ratio a 
for gene networks with 60 genes with k/N = 0.05. Both input and output data are subject to 
observational noise of r\ = 0.30. The performance was obtained using both Wu (open square) 
and Wl2 (filled circle) prescriptions and a = 0.01. Top: the a priori probabilities are equal, 
i.e. pfj = p^j = = 1/3 for all i and j. Medium: the a priori probabilities are set to be 
pfj = p^j = 0.025 and p?- = 0.95 for all i and j. Bottom: the a priori probabilities are set to be 
Pin = Pij = 0-01 an d Pi~ = 0.98 for all i and j. 
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FIG. 4: (Color online) Subnetwork identification. A: Information gain I g obtained for three "over- 
expression experiments". Firstly, the gene which belongs to the subnetwork of interest, gene gl, 
is initially over expressed, then the input-output network is measured, this measurement is subject 
to observational noise of rj = 0.30. The information gain of this experiment is computed for each 
gene and the genes with I g greater than a given threshold, are selected. I g suggests that gene g6 is 
regulated by gl. Repeating the above step with gene g6, it appears that genes g2, g3 are regulated 
by g6. The above step is repeated with gene g2 and subsequent genes with high information gain 
values in subsequent experiments. B: List of experiments, the first column corresponds to the gene 
which was overexpressed in each experiment, the second column corresponds to the genes which 
appear to be regulated by the overexpressed gene. C: Subnetwork inferred 23 out of 24 interactions 
correctly (solid edges) by this inference procedure using Wl2 prescription and the above ten "over- 
expression experiments" together. The information contained in list B was included as a priori 
probabilities, i.e., they were set pfj = p~j = 0.5 and p\j = 0.0 for all i,j pairs indicated in the list, 
and pf; = p~, = 0.01 otherwise. 
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